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1 .  Introduction 

One  of  the  major  problems  that  arises  in  the  statistical  analysis  of 
simulation  output  is  the  generation  of  confidence  intervals  for  parameters 
of  interest.  However,  a  major  practical  obstacle  remains:  Coverage  rates 
tend  to  be  substantially  lower  than  the  confidence  level  indicated.  This 
phenomenon  manifests  itself  even  in  those  cases  in  which  the  procedures 
have  an  asymptotically  consistent  large  sample  theory.  For  a  discussion  of 
this  problem,  we  refer  the  reader  to  Section  8.5.1  of  Law  and  Kelton 
(1982). 

In  this  paper,  we  will  study  asymptotic  expansions  associated  with  the 
error  of  the  coverage  probability.  We  begin,  in  Section  2,  with  a 
discussion  relating  to  the  interpretation  of  a  confidence  interval.  The 
asymptotics  associated  with  confidence  interval  error  turn  out  to  be 
interpretation-dependent.  In  Section  3,  we  consider  error  expansions  for 
confidence  intervals  for  the  mean  of  a  sequence  of  independent  and 
identically  distributed  (i.i.d.)  random  variables.  This  case  is  of 
interest  when  the  simulator  is  concerned  with  analyzing  the  output  of  a 
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terminating  simulation.  It  is  also  intimately  connected  to  the  application 
of  the  technique  of  replication  to  the  steady-state  confidence  interval 
problem.  In  this  context,  we  are  able  to  precisely  Identify  the  primary 

t 

sources  of  error  in  the  converage  rate,  and  we  discuss  the  relevance  of  the 
error  expansion  to  the  choice  between  using  a  t-variate  or  normal  in  the 
confidence  interval  procedure.  In  Section  4,  we  examine  the  regenerative 
confidence  interval  that  arises  when  a  simulator  is  interested  in  the 
steady-state  mean  of  a  regenerative  stochastic  process.  We  show  that  the 
coverage  error  has  a  form  similar  to  that  in  the  i.i.d.  case.  Our  results 
VMuicate  that  the  coverage  rate  difficulty  inherent  in  the  regenerative 
confidence  interval  is  of  the  same  order  of  magnitude  as  the  coverage  rate 
error  faced  in  the  classical  i.i.d.  case.  Finally  in  Section  5,  we  state 
some  conclusions,  and  briefly  discuss  two  new  small  sample  confidence 
interval  methodologies  that  appear  to  have  enhanced  coverage  properties. 


2.  The  Interpretation  of  a  Confidence  Interval 

Let  {Xn;  n  >_  0}  be  a  sequence  of  random  variables  representing  the 
output  of  a  simulation,  and  suppose  that  a  confidence  interval  for  the 
parameter  0  is  required.  The  goal  of  the  simulator  is  to  find  a  sequence 


of  random  variables 


L  =  L  (Xn,  ...,  X  ) 
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such  that  the  random  interval  [Ln,Rn]  corresponds  to  a  "confidence" 
set  for  9*  In  precise  terms,  [Ln,RnJ  is  a  100(1-a)*  confidence 
interval  for  9  if 

(2.1)  P{6  e  CLn,RnD>  =  1-a  . 

However,  (2.1)  does  not  fully  specify  the  choice  of  and  Rn,  as 
the  following  example  shows. 

'?.?)  Example.  Let  {Xn;  n  )  0}  be  normal  random  variables  wxth 
unknown  mean  6  and  known  variance  o^,  denoted  N(6,o^).  Let 
®(x)  =  P{N(0,1)  £  x},  and  put  z(p)  =  ®“^(p)  for  0  <  p  <  1.  Then,  it 
is  easily  verified  that 

[Xn  -  z(p+1-a)o/n^2,  Xn  -  z(p)c/n^2]  , 

where  X^/n,  is  a  100(1-a)%  confidence  interval  for  9,  provided 

0  <  p  <  a. 

In  the  above  example,  most  simulators  would  agree  that  the  "correct" 
choice  of  p  is  a/2.  There  are  two  reasons  for  that  choice.  First,  the 
length  of  the  interval  is  minimal  for  p  =  a/2.  Secondly,  the  interval 
[Ln,Rn]  obtained  through  that  choice  has  the  property  that 

P(0  >  L  }  =  P{9  <  R  }  =  1-a/2  , 

—  n  —  n 

P{9  <  L  >  *  P{9  >  R  >  =  a/2  . 
n  n 


(2.3) 
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We  shall  call  an  interval  [Ln,Rn]  satisfying  (2.3)  a  100(1-ct)S 
balanced  confidence  interval  for  8.  Observe  that  (2.3)  implies  (2.1)  so 
that  any  balanced  confidence  interval  is  a  confidence  interval  in  the  sense 
of  (2.1). 

In  many  simulation  problems,  it  appears  that  the  simulation  practi¬ 
tioner  would  have  a  preference  for  a  balanced  interval.  For  example, 
suppose  that  a  simulation  of  a  queueing  system  produces  a  90S  confidence 
interval  [Ln,Rn]  for  the  mean  customer  waiting  time  6.  The  simu¬ 
lator  can  conclude  only  that 

0  <  P{0  <  Rn),  P{0  >  Ln)  <  0.1 
whereas,  for  a  balanced  interval,  the  conclusion 

P(8  >  R „}  =  P{9  <  L  }  =  0.05 
n  n 

is  possible.  Clearly,  a  balanced  interval  gives  more  information  to  the 
simulator. 

It  should  be  pointed  out  that  much  of  the  statistical  theory  of  confi¬ 
dence  intervals  ln  a  parametric  setting  relates  to  balanced  intervals.  For 
a  discussion  of  the  desirability  of  balanced  intervals  from  a  Bayesian 
viewpoint,  see  Efron  (1981). 
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3.  Coverage  Rates  for  I.I.P.  Confidence  Intervals 

Let  {Xn;  n  0}  be  a  sequence  of  i.i.d.  random  variables  with 
0  <  =  var(Xn)  <  Let  us  first  examine  the  results  in  the  simplest 

possible  setting:  The  goal  is  to  produce  a  confidence  interval  for 
p.  =  EXn,  given  that  is  known. 

The  standard  confidence  interval  procedure,  in  this  case,  starts  from 
a  Central  Limit  Theorem  (CLT)  for  n^/2  (Xn-p)/a.  Simple  algebraic 
manipulation  shows  that  [Ln(p),  R^p)]  is  an  approximate  100(1-a)X 
confidence  interval  for  p,  where 


Ln(p)  =  Xn  -  z(p+1-a)a/n1/2 
Rn(p)  =  Xn  -  z(p)o/n1^2 


for  0  <  p  <  a. 

oefore  proceeding,  let  us  observe  that  the  coverage  rate  error  in  the 
interval  [*-n(p)»  Rn(p)],  denoted  en(p)»  is  given  by  jP{p  e  [Ln(p),Rn(p)]> 
-  (1-o)J.  As  for  the  balanced  situation,  it  is  clear  that  the  interval 
[Ln(p),  Rn(p)]  is  asymptotically  balanced  only  if  p  =  a/2.  Thus,  we 
shall  henceforth  restrict  our  discussion  of  balanced  error  to  this  case, 
and  designate  (en,en>  as  our  error  descriptor,  where 
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e*  *  |P{p  <  Ln(a/2)>  -  «/2| 
ejj  -  (P<ii  >  Rn(«/2)}  -  a/2 1 

Our  error  estimates  will  be  written  in  terms  of  the  "big  0",  "little  o” 
notation  (i.e.,  g(n)  =  0(f(n))  if  g(n)/f(n)  is  bounded;  g(n)  =  o(f(n)) 
if  g (n)/f(n)  goes  to  zero). 

&  I  M 

(3.1)  Theorem,  (i)  If  EX„  <  •  ,  then  the  error  terms  e  (p),  e  ,  e  are 

n  n  n  n 

all  0(n-!/2). 


(ii)  If,  in  addition,  Xn  has  a  distribution  with  a  Lebesgue 
density  component,  then 

a)  en(p)  a  JSk(XQ)-(g(z(p))  -  g(z(p+1>a))) |/6n^2  +  0(-^)  , 

b)  e*  =  |Sk(X0).g(z(a/2))|/6n1/2  +  0(-J)  =  e*  , 

where 

g(x)  *  (1-x2) d/dx(*(x))  ,  and  Sk(XQ)  a  E(X0-p)3/o(X0)3  . 


Proof.  Part  (i)  of  the  theorem  is  a  direct  consequence  of  the  Berry-Esseen 
theorem  (see  Feller  (1971),  p.  542).  The  proof  of  the  second  half  of  the 
result  follows  from  simple  algebraic  manipulation  of  the  expansions  cited 
in  Theorems  1  and  3  of  Section  16.4  of  Feller  (1971). 


The  first  part  of  the  theorm  tells  us  that  under  an  appropriate  moment 

•  f* 

assumption',  the  errors  en(p) ,  en,  and  en  decrease  at  least  as  fast  as 
n~1/2,  for  all  p.  The  second  half  of  the  result  shows  that  when  Xn 

O 

has  a  density  component,  then  en(a/2)  is  0(1/n),  whereas  for  p  *  a/2, 

- 1/2 

en(p)  is  0(n  ).  On  the  other  hand,  the  balanced  error  for  [Ln(a/2), 

Rn(a/2)]  Is  always  0(n~^2). 

Thus,  using  p  =  a/2  in  the  interval  [Ln(p),  R„(p)]  produces  a 

confidence  interval  with  the  "best"  coverage  rate  (i.e.,  0(1/n)  as  opposed 

to  0(n“^2)),  as  well  as  the  asymptotically  minimal  length.  However, 

c^a/2),  Rh ( a/2)  J  achieves  this  result  in  a  very  curious  way.  Part 

(ii.b)  of  Theorem  3.1  shows  that  P{p  £  ^(0/2)}  will  differ  from  a/2 
-1/2 

by  0(n  ).  Hence,  if  P{p  £  Ln(a/2)}  is  greater  (say)  than  all  by 

0(n”^),  then  P{p  _>  Rn(a/2)}  must  be  less  than  a/2  by  precisely  the 
same  0(n-1/2)  error  (in  order  that  en(a/2)  be  0(1/n)).  This 
suggests,  in  this  case,  that  a  "better"  confidence  interval  would  be 
achieved  by  shifting  the  interval  to  the  left  slightly.  It  Is  also  worth 
observing  that  the  coverage  rate  accuracy  of  p  s  a/2  is  highly  unstable 
in  the  sense  that  en(a/2+q)/en(a/2)  «  for  any  non -zero  t). 

Some  caution  should  be  exercised  in  trying  to  extend  the  second 
conclusion  of  the  theorem  to  the  case  of  discrete  random  variables.  The 


result  depends  on  an  Edgeworth  expansion  which  breaks  down  in  the  discrete 
case  (see  [3],  p.  539,  for  a  related  comment). 

Our  above  analysis  required  that  a2  be  known.  Of  course,  in 
general  the  simulator  does  not  know  this  parameter,  and  hence  it  must  be 
estimated  from  the  simulation  output  sequence  {Xn}.  The  key  result  in 


forming  confidence  intervals  in  this  setting  is  the  CLT  for  tn,  where 
tn  =  n^2(Kn-|i)/sn 

s2  1  l  Xf  -  (X  )2  . 

n  n-1  i  n 

The  intervals  [Ln(p),  Rn(p)]  are  approximate  100(1-a)S  confidence 
intervals  for  0  <  p  <  a,  provided  Ln(p)  =  Xn-z(p+1-a)sn/n^2, 

Rn(p)  =  Vz(p)sn/f,1/2* 

8  9  p 

v3 .£)  Theorem,  (i)  If  EX^  <  •»,  then  the  error  terms  en(p)»  En»  en  are 
all  (Kn-1^). 

(ii)  If,  in  addition,  Xn  has  a  distribution  with  a  Lebesgue 
density  component,  then 

a)  en(p)  =  |slc^xo^ * <h(z(p))  -  h(z(p+1-a)))  |/n^2  +  0(1)  , 

b)  e*  =  |Sk(X0).h(z(a/2))/n1/Z  ♦  0(1)  =  , 

where  h(x)  =  [(1+2x)d/dx(®(x))]/6. 

The  proof  of  this  theorem  is  similar  to  that  of  Theorem  3.1.  The  key 
step  is  to  prove  a  Berry-Essen  type  result  and  obtain  an  Edgeworth 
expansion  for  tn;  this  can  be  found  in  Glynn  (1982). 


Notice  that  the  conclusions  of  this  theorem  are  qualititively  very 

similar  to' those  of  Theorem  3.1,  where  a 2  was  known.  For  example, 

tla/2)  is  again  0{1/n),  whereas  ef,  er  are  both  0(n’^2).  However, 
n  n  n 

the  coefficient  in  n"1/2  has  changed.  The  proof  of  the  Edgeworth 

expansion  mentioned  above  shows  that  the  error  of  order  n"1/2  occurs 

precisely  as  a  consequence  of  the  skewness  Sk(Xg)  of  {X^ ;  n  ^  0),  and 

the  correlation  between  s2  and  X  .  This  is  in  contrast  to  the  situation 

n  n 

where  a 2  is  known,  in  which  case  all  the  n“^/2  error  emanates  from 

skewness  alone.  Recall  that  the  coefficient  of  skewness  Sk(Xg) 

..icoiures  the  asymmetry  of  a  distribution. 

Before  concluding  this  section,  we  turn  to  the  question  of  using 

t-variates  rather  than  a  normal  distribution  to  generate  confidence 

intervals.  Let  T^  have  a  Student's  t-distribution  with  k  degrees  of 

freedom.  Put  rk<x)  =  P{Tk£x},  and  set  zk(p)  =  rj^ 1  ( p ) -  Many  simulators 

(e.g.,  [6],  p.  288)  advise  the  use  of  the  interval  [L^fp),  fin(p)],  where 

L  (p)  =  X-z  .(p+1-a)s  /n^2,R  (p)  =  X-z  ,(p)s  /n1^2,  rather  than  the 
n  n  n-i  n  n  n  n- l  r  n 

*  -JL  *r 

previously  defined  interval  for  p.  Let  e  (p),  c  ,  e  be  the  coveraqe 

n  n  n 

rate  error  and  balancing  errors  associated  with  the  t-variate  confidence 


interval. 
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(3.3)  Theorem.  For  0  <  p  <  a, 

(l)  en(p)  3  en^  *  °^n_1^2^ 

(li)  ij(p)  =  e*(p)  ♦  o(n‘1/2) 

ejj(p)  =  e£(p)  ♦  o(n"^2) 

The  proof  of  this  result  can  be  found  in  [4].  The  content  of  this 
i-h^rem  is  the  that  t-variate  modification  is  of  "small  order"  in  the  sense 
that  the  leading  term  in  the  error  (of  order  n“^)  is  precisely  the 
same  as  that  obtained  via  normal  theory.  However,  it  should  be  mentioned 
that  the  coverage  error  en(a/2)  may  differ  from  en(a/2)  in  its  leading 
term.  In  any  case,  it  seems  clear  that  the  desirability  of  using  a 
t-variate,  as  opposed  to  a  normal,  deserves  further  study. 


4.  Coverage  Rates  for  Regenerative  Confidence  Intervals 

Let  (X^;  t  >_  0)  be  a  regenerative  stochastic  process.  Then,  there 
exist  random  times  T-|,  Tj,  ...  such  that  for  any  (measurable) 
real-valued  function  f,  {(Y^,  k  ^  1}  is  i.i.d.,  where 


k+1 

\  =  /  f(X,)  ds 

% 


tk  =  Tk+1  '  Tk  * 


1 

i 
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The  goal  of  the  simulator  is  to  find  a  confidence  interval  for  the 
steady-state  mean  of  { f ( ) ;  t  ^  0).  It  can  be  shown  that  this  is 
equivalent  to  obtaining  an  interval  for  r  =  EYi/Evj  (see  Crane  and 
Iglehart  (1975)).  The  regenerative  method  for  output  analysis  depends  on  a 
CLT,  under  the  assumption  0  <  a^Y^-rx^)  <  ®,  Ex^  <  ®,  for  the  statistic 
n1/2(rn_r)/Vnj  Where 


n 

rn  =  l 
n  kul 


V  I 

k=1 


V  = 


<VT1> 


n 

(  l 

k=1 


Y,  -  2r 


n 

n  k=1 


?  n 

YkTk  +  rn  l 
K  K  n  k=1 


As  for  the  i.i.d.  case  of  Section  3,  the  normal  approximation  yields 
[Ln(p) ,  Rn(p)]  as  an  approximate  100(1~a)%  confidence  interval  for 
r,  where  L^p)  =  rn-z(p+1-a)vn/n^2,  R^p)  =  r^zfpW^/n1^2. 

8  8 

(4.1)  Theorem,  i)  If  EY  <  ®  and  Ex  <  ®,  then  the  error  terms 

n  n 

en(p),  €*,  are  all  0(n-1^2). 


(ii)  If,  in  addition,  the  distribution  of  (Yn,xn)  has  a 
component  which  is  a  Lebesgue  density  in  the  plane,  then 


a)  en(p)  =  Jk(z(p) )  -  k(z(p+1-a) )  |  /n^2  +  0(1/n). 

b)  ej  =  |k(z(a/2))|  ♦  0(1/n)  = 


where 
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k(x)  =  (a+px2)  d/dx(«(x))  , 
a  =  E 1^/60*  {1 J 

0  =  EZ^c3^)  -  Et1Z1/cj(Z1)Et1 


This  theorem  is  a  consequence  of  a  Berry -Esseen  result  and  Edgeworth 
expansion  for  regenerative  processes,  which  appears  in  [4]. 

The  important  point  here  is  the  similarity  to  the  i.i.d.  case.  The 
order  of  error  is  precisely  the  same  as  that  which  we  encountered  in 
Section  3.  Hence,  it  could  be  argued,  on  the  basis  of  error  comparison, 
that  the  steady-state  simulation  problem  is  no  more  difficult  than  the 
terminating  simulation  problem,  provided  that  the  regenerative  method  is 
used. 


As  for  the  form  of  the  error  coefficient  in  we  observe  that 

Sk(Zi)  plays  an  important  role,  as  in  the  i.i.d.  case.  Additional  terms 
in  Et^Z^/Et^ »c(Z^ )  also  appear,  however:  These  are  contributed  by  the 
bias  of  rn.  Thus,  the.  prime  sources  of  error  in  the  regenerative  case 
are  asymmetry  effects,  correlation  between  point  and  variance  estimates, 
and  bias  problems  in  the  point  estimate.  This  has  an  important  implication 
for  research  efforts  directed  at  producing  "correct"  coverage  rates  for 
confidence  intervals.  It  shows  that  reducing  one  source  of  error,  such  as 
bias  in  the  point  estimate,  should  not  be  expected  to  necessarily  improve 


coverage 
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5.  Concluding  Remarks 

In  this  paper,  we  have  provided  an  overview  of  the  coverage  error  prob¬ 
lem  for  confidence  interval  generation  in  simulatxo>*.  We  have  shown  that 
the  qualitative  character  of  the  error  appears  reasonably  insensitive  to 
the  "fine"  structure  of  the  simulation  output  sequence.  The  regenerative 
case  indicates  that  point  estimate  bias,  correlation  between  point  and 
variance  estimators,  and  asymmetry  all  play  equally  important  roles  in 
determining  coverage  rates  in  the  steady-state  simulation  problem.  Our 
results  also  show  the  need  for  more  research  on  the  question  of  whether 
i -variates  or  normals  should  be  used  to  generate  intervals.  Furthermore, 

.ve  have  shown  that  the  parameter  p  plays  a  critical  role  in  determining 
the  amount  of  error.  However,  we  caution  the  reader  that  although  p  =  a/2 
appears  optimal  in  the  sense  of  error,  the  result  is  highly  unstable. 
Finally,  we  have  introduced  the  concept  of  a  balanced  confidence  interval 
and  shown  that  its  error  asymptotics  are  somewhat  different  from  the 
standard  error  criterion. 

In  [4],  we  study  two  procedures  that  appear  to  have  promising  coverage 
characteristics.  The  first  technique  is  a  regenerative  bootstrap  (see  [2] 
for  the  bootstrap  in  the  i.i.d.  case),  and  the  second  method  involves  an 
application  of  a  so-called  "normalization"  transformation.  The  latter 
procedure  is  based  on  an  idea  of  Oohnson  (1978).  The  asymptotic  error 
expansion  for  these  intervals  indicates  an  improvement  over  currently  used 
intervals.  These  improvements  have  also  manifested  themselves  empirically. 
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